//#include <vld.h>
#include <fstream>
#include <iostream>
#include <deque>

//#include "cut.h"
#include "augmentedctree.h"

using namespace PBVP;

AugmentedContourTree::AugmentedContourTree(ContourTree* conTree)
{
	m_pIntvTree = NULL;
	m_pArcs = new SuperArcMap;
	augment(conTree);
	reduce();
	finalize();
}

AugmentedContourTree::AugmentedContourTree()
{
	m_pIntvTree = NULL;
	m_pArcs = new SuperArcMap;
}

AugmentedContourTree::~AugmentedContourTree()
{
	if(m_pIntvTree != NULL) delete m_pIntvTree;
	delete m_pArcs;
}

void AugmentedContourTree::augment(ContourTree* conTree)
{
	int i, j;
	int nv = (int)conTree->numOfNodes();
	int *x = new int[nv];
	int *b = new int[nv];
	int *p_map = new int[nv];
	for (i = 0; i < nv; i++)
	{
		x[i] = 0;
		b[i] = 0;
	}

	std::deque<int> leafs;
	// iterate through the contour tree
	i = 0;
	std::vector<CriticalPoint>* pCPVector = conTree->getSortedPoints();
	std::vector<CriticalPoint>::iterator cit = pCPVector->begin();
	while(cit != pCPVector->end())
	{
		SuperNode& node = (*cit);
		int id = node.id;
		p_map[id] = i;
		this->addNode(node);
		if(conTree->isLeaf(id))
		{
			leafs.push_back(id);
		}
		i++;
		++cit;
	}

	while (!leafs.empty())
	{
		int nid;
		int id = leafs.front();
		leafs.pop_front();
		SuperNode& node = conTree->getNode(id);
        int degree = node.degree();
		if (degree > 0)
		{ // not the last node
			nid = node.getNeighborAt(0);
			int x1 = node.getTagAt(0);
			SuperNode& tmp = conTree->getNode(nid);
			int x2 = tmp.getNeighborTag(id);
			i = p_map[id]; 
			j = p_map[nid];

			int LSi = node.LS;
			int USi = node.US;
			int dbei = node.dbe;

			int delta = (i < j)? 1:-1;
			int xe = delta * (x[i] - USi + LSi);
			int be = delta * (b[i] + dbei);
			x[j] += delta * xe;
			b[j] += delta * be;
#ifdef _DEBUG
//printf("add acr (%d %d) with xe = %d, be = %d\n", id, nid, xe, be);
#endif
			addArc(id, nid, xe, be, -1, x1, x2);
		}
		conTree->removeNode(id);
		if (degree > 0)
		{
			if (conTree->isLeaf(nid)) leafs.push_back(nid);
		}
	}
	delete[] x;
	delete[] b;
	delete[] p_map;
	delete pCPVector;

	setMaxID(conTree->getMaxID());
}

void AugmentedContourTree::calcBettiNumbers()
{
	SuperArcMap::iterator it = m_pArcs->begin();
	while(it != m_pArcs->end())
	{
		SuperArc& arc = (*it).second;
		arc.calcBettiNumber();
		++it;
	}
}

void AugmentedContourTree::reduce()
{
	printf("numbe of actree nodes before reduce: %d\n", numOfNodes());
	SuperNodeMap::iterator it = nodes.begin();
	while (it != nodes.end())
	{
		SuperNode &node = (*it).second;
		int id = node.id;
        int LS = node.LS;
		int US = node.US;
		int dbe = node.dbe;
		++it;
		if (node.degree() == 2 && LS == 0 && US == 0)
		{
			removeNode(id);
		}
		// printf("number of nodes: %d\n", nodes.size());
	}
	printf("number of actree nodes after reduce: %d\n", numOfNodes());
}

void AugmentedContourTree::removeNode(int n)
{
	SuperNode& node = nodes[n];

	assert(node.degree() <= 2);	// only node with degree 1 or 2 can be removed
	if (node.degree() == 1)
	{
		int nid = node.getNeighborAt(0);
		SuperNode& nei = nodes[nid];
		nei.removeNeighbor(n);
		nei.removeArc(n);
		int arcid = node.getArc(nid).id;
		m_pArcs->erase(arcid);
	} else if (node.degree() == 2)
	{
		int j = node.getNeighborAt(0);
		int k = node.getNeighborAt(1);
		int xe, be;
		SuperArc& arc1 = node.getArc(j);
		SuperArc& arc2 = node.getArc(k);
		if (arc1.betti1() != arc2.betti1() || arc1.betti2() != arc2.betti2())
		{
			// should not delete the node
			return;
		}
		xe = arc1.xe;
		be = arc1.be;

		SuperNode& nei1 = nodes[j];
		SuperNode& nei2 = nodes[k];
		
		m_pArcs->erase(arc1.id);
		m_pArcs->erase(arc2.id);
		nei1.removeNeighbor(n);
		nei1.removeArc(n);
		nei2.removeNeighbor(n);
		nei2.removeArc(n);

		// add a new arc between the two neighbors	
		addArc(j, k, xe, be, n);
	}
	nodes.erase(n);
}

void AugmentedContourTree::addArc(int v1, int v2, int xe, int be, int nid, int x1, int x2)
{
	SuperNode &node1 = nodes[v1];
	SuperNode &node2 = nodes[v2];

	node1.addNeighbor(v2, x1);
	node2.addNeighbor(v1, x2);

	SuperArc sa(v1, v2, xe, be);
	if (nid >= 0) sa.addIntraNode(nid);

	(*m_pArcs)[sa.id] = sa;
	node1.addArc(v2, sa);
	node2.addArc(v1, sa);
}

bool AugmentedContourTree::saveFile(const char* fname)
{
	//std::vector<SuperNode>* sortedNodes = getSortedPoints();
	std::ofstream out(fname);
	out << 1 << std::endl;
	int nc = (int)nodes.size(), na = (int)(m_pArcs->size());
	out << nc << " " << na << std::endl;

	SuperNodeMap::iterator cit = nodes.begin();
	while ( cit != nodes.end() )
	//std::vector<SuperNode>::iterator cit = sortedNodes->begin();
	//while( cit != sortedNodes->end() )
	{
		SuperNode& node = (*cit).second;
		//SuperNode& node = nodes[(*cit).id];
		out << node.id << " " << node.val << " ";
		out << node.degree() << " ";

		for (int i = 0; i < node.degree(); i++)
		{
			out << node.getNeighborAt(i) << " " << node.getTagAt(i) << " ";
		}
		out << std::endl;
		++cit;
	}
	//delete sortedNodes;

	// output arcs
	SuperArcMap::iterator ait = m_pArcs->begin();
	while (ait != m_pArcs->end() )
	{
		SuperArc &arc = (*ait).second;
		out << arc.v1 << " " << arc.v2 << " ";
		out << arc.xe << " " << arc.be << "  ";
		out << 1 << " " << arc.intra_node << std::endl;
		++ait;
	}
	out.close();
	return true;
}

bool AugmentedContourTree::loadFile(const char* fname)
{
	int nc, na, i;

	std::ifstream in(fname);
	if (!in.is_open())
	{
		std::cerr << "cannot open file " << fname << std::endl;
		return false;
	}
	int head;
	in >> head;
	if (head != 1)
	{
		std::cerr << "not an augmented contour tree, aborted\n";
		in.close();
		return false;
	}
	in >> nc; in >> na;
	for (i = 0; i < nc; i++)
	{
		int id, nid, xid, degree;
		float val;
		in >> id;
		in >> val;
		in >> degree;
		addNode(id, val);
		for (int j = 0; j < degree; j++)
		{
			in >> nid;
			in >> xid;
			nodes[id].addNeighbor(nid, xid);
		}
	}

	for (i = 0; i < na; i++)
	{
		int v1, v2, xe, be, count, intra_node;
		in >> v1;
		in >> v2; 
		in >> xe;
		in >> be;
		in >> count;
		in >> intra_node;
		addArc(v1, v2, xe, be, intra_node);
	}
	in.close();
	finalize();
	return true;
}

void AugmentedContourTree::finalize()
{
	int i, j; 
	int	nv = (int)nodes.size();
	float *vals = new float[nv];

	printf("construct interval tree...\n");
	// traverse the critical point list
	std::vector<CriticalPoint> *cpList = getSortedPoints();
	std::vector<CriticalPoint>::iterator cit = cpList->begin();
	i = 0;
	while ( cit != cpList->end() )
	{
		vals[i] = (*cit).val;
		++cit;
		i++;
	}
	// compacting critical vals
	for (i = 1, j = 0; i < nv; i++)
	{
		if (vals[i] != vals[j])
		{
			vals[++j] = vals[i];
		}
	}

	if(m_pIntvTree != NULL) delete m_pIntvTree;
	m_pIntvTree = new IntervalTree<SuperArc>(j+1, vals);
	delete[] vals;
	delete cpList;
	
	// add arcs to the interval tree
	SuperArcMap::iterator ait = m_pArcs->begin();
	// LL_Node<SuperArc>* arc_ptr = arcs.head();
	while( ait != m_pArcs->end() ) {
		int n1 = (*ait).second.v1;
		int n2 = (*ait).second.v2;
		float x1 = nodes[n1].val;
		float x2 = nodes[n2].val;
		m_pIntvTree->insert((*ait).second, std::min(x1, x2), std::max(x1, x2));
		++ait;
	}
	m_pIntvTree->done();
/*	commented by xiaoyu 3/25/08
	// set IDs the super arcs 
	i = 0;
	// arc_ptr = arcs.head();
	QList<SuperArc>::iterator ait = arcs.begin();
	while ( ait != arcs.end() )
	{
		(*ait).id = i;
		i++;
		++ait;
	}  */
}


BettiNumber AugmentedContourTree::getBettiNumbers(float x)
{
	BettiNumber betti;

	printf("get Betti number of %f\n", x);
	if(m_pIntvTree == NULL) {
		std::cerr << "construct interval tree before querying Betti Numbers\n";
		return betti;
	}
	vector< SuperArc >* p_results = m_pIntvTree->query(x);
	vector< SuperArc >::iterator it;
	int count = 0;
	for(it = p_results->begin(); it != p_results->end(); ++it) {
		SuperArc& arc = (*it);
		betti.b0 += 1;
		betti.b1 += arc.betti1();
		betti.b2 += arc.betti2();
		count ++;
		printf("%d th arc %d: (%f, %f) \n", count, arc.id, nodes[arc.v1].val, nodes[arc.v2].val);
	}
	delete p_results;
	return betti;
}


std::vector<SuperArc>* AugmentedContourTree::cut(float x)
{
	if(m_pIntvTree == NULL) {
		std::cerr << "need to construct the interval tree first\n";
		return NULL;
	}

	std::vector<SuperArc>* pArcs = m_pIntvTree->query(x);	
	return pArcs;
}
